Objectives

“Identify a single vs. a multi-band raster file.”
“Import multi-band rasters into R using the raster package.”
“Plot multi-band color image rasters in R using the ggplot package.”

Keypoints

“A single raster file can contain multiple bands or layers.”
“Use the stack() function to load all bands in a multi-layer raster file into R.”
“Individual bands within a stack can be accessed, analyzed, and visualized using the same functions as single bands.”

If you’re just opening up this project, you’ll have to load the packages and data used in previous episodes.

# load packages 
library(raster)
library(rgdal)
library(ggplot2)
library(dplyr)

Getting Started with Multi-Band Data

In this episode, the multi-band data that we are working with is imagery collected using the NEON Airborne Observation Platform high resolution camera over the NEON Harvard Forest field site. Each RGB image is a 3-band raster. The same steps would apply to working with a multi-spectral image with 4 or more bands - like Landsat imagery.

If we read a RasterStack object into R using the raster() function, it only reads in the first band.

# import RGB using raster function
# only imports first band
rgb_band1 <- raster("data/raster/HARV_RGB_Ortho.tif")

# look at some metadata
rgb_band1

## class      : RasterLayer 
## band       : 1  (of  3  bands)
## dimensions : 1571, 1938, 3044598  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho 
## values     : 0, 255  (min, max)

We can use the raster() function to import specific bands in our raster object by specifying which band we want with band = N (N represents the band number we want to work with). To import the green band, we would use band = 2.

# import using raster, specify band number
rgb_band2 <-  raster("data/raster/HARV_RGB_Ortho.tif", band = 2)

# look at some metadata
rgb_band2

## class      : RasterLayer 
## band       : 2  (of  3  bands)
## dimensions : 1571, 1938, 3044598  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho 
## values     : 0, 255  (min, max)

Raster Stacks

Next, we will work with all three image bands (red, green and blue) as an R RasterStack object. We will then plot a 3-band composite, or full color, image.

To bring in all bands of a multi-band raster, we use thestack() function.

# import all bands using stack
rgb_stack <- stack("data/raster/HARV_RGB_Ortho.tif")

# look at some metadata
rgb_stack

## class      : RasterStack 
## dimensions : 1571, 1938, 3044598, 3  (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## names      : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3 
## min values :                0,                0,                0 
## max values :              255,              255,              255

We can view the attributes of each band in the stack in a single output. Or, if we have hundreds of bands, we can specify which band we’d like to view attributes for using an index value:

# look at information about the bands
rgb_stack@layers

## [[1]]
## class      : RasterLayer 
## band       : 1  (of  3  bands)
## dimensions : 1571, 1938, 3044598  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.1 
## values     : 0, 255  (min, max)
## 
## 
## [[2]]
## class      : RasterLayer 
## band       : 2  (of  3  bands)
## dimensions : 1571, 1938, 3044598  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.2 
## values     : 0, 255  (min, max)
## 
## 
## [[3]]
## class      : RasterLayer 
## band       : 3  (of  3  bands)
## dimensions : 1571, 1938, 3044598  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.3 
## values     : 0, 255  (min, max)

rgb_stack[[2]]

## class      : RasterLayer 
## band       : 2  (of  3  bands)
## dimensions : 1571, 1938, 3044598  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732483, 4712956, 4713349  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.2 
## values     : 0, 255  (min, max)

We can also use the ggplot functions to plot the data in any layer of our RasterStack object. Remember, we need to convert to a data frame first.

Each band in our RasterStack gets its own column in the data frame. Note the column names are from the name of the file with the band number.

# convert to data frame for plotting
rgb_stack_df  <- as.data.frame(rgb_stack, xy = TRUE)

# look at column names
str(rgb_stack_df)

## 'data.frame':    3044598 obs. of  5 variables:
##  $ x               : num  731999 731999 731999 731999 732000 ...
##  $ y               : num  4713349 4713349 4713349 4713349 4713349 ...
##  $ HARV_RGB_Ortho.1: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ HARV_RGB_Ortho.2: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ HARV_RGB_Ortho.3: num  0 0 0 0 0 0 0 0 0 0 ...

Let’s create a histogram of the first band and a raster plot of the second band.

# histogram of band 1
ggplot() +
  geom_histogram(data = rgb_stack_df, aes(HARV_RGB_Ortho.1))

# plot band 2
ggplot() +
  geom_raster(data = rgb_stack_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho.2)) + 
  coord_quickmap()

Create A Three Band Image

To render a final three band, colored image in R, we use the plotRGB() function.

This function allows us to:

Identify what bands we want to render in the red, green and blue regions. The plotRGB() function defaults to a 1=red, 2=green, and 3=blue band order. However, you can define what bands you’d like to plot manually. Manual definition of bands is useful if you have, for example a near-infrared band and want to create a color infrared image.
Adjust the stretch of the image to increase or decrease contrast.

Let’s plot our 3-band image. Note that we can use the plotRGB() function directly with our RasterStack object (we don’t need a dataframe as this function isn’t part of the ggplot2 package).

# plot RGB image
plotRGB(rgb_stack,
        r = 1, g = 2, b = 3)

The image above looks pretty good. We can explore whether applying a stretch to the image might improve clarity and contrast using stretch="lin" or stretch="hist".

# edit plot using a linear stretch
plotRGB(rgb_stack,
        r = 1, g = 2, b = 3,
        scale = 800,
        stretch = "lin")

The material has been adapted from The Carpentries Introduction to Geospatial Raster and Vector Data with R lesson licensed under CC-BY 4.0

Work with multi-band rasters

Getting Started with Multi-Band Data

Raster Stacks

Create A Three Band Image